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Abstract. We present a Markov chain Monte Carlo technique for detecting 
gravitational radiation from a neutron star in laser interferometer data. The 
algorithm can estimate up to six unknown parameters of the target, including 
the rotation frequency and frequency derivative, using reparametrization, delayed 
rejection and simulated annealing. We highlight how a simple extension of the 
method, distributed over multiple computer processors, will allow for a search 
over a narrow frequency band. The ultimate goal of this research is to search 
for sources at a known locations, but uncertain spin parameters, such as may be 
found in SN1987A. 



1. Introduction 

Rapidly rotating neutron stars could be an important source of gravitational wave 
signals. Several mechanisms have been proposed that would cause them to emit quasi- 
periodic gravitational waves 

Interferometric gravitational wave detectors are now operating in numerous 
locations around the world ^JEHIHE], and much work has gone into the development of 
dedicated search algorithms for these signals. Radio observations can provide the sky 
location, rotation frequency and spin-down rate of known pulsars, and this knowledge 
simplifies the analysis. This was the case for the recent search for a signal from 
PSR J1939+2134 T\. When the position and phase evolution of a source are not 
known, all-sky hierarchical strategies are required, and these have huge computational 
requirements [HHH]. 

Here we concentrate on the search for a gravitational wave signal from a known 
location, but where spin parameters of the rotating neutron star are not well known 
(but within a narrow band). SN1987A is a good example of a poorly parameterised 
source for which the sky location is known, but where there are large uncertainties in 
the frequency and spin-down parameters of a putative neutron star ^Hj- In particular, 
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we consider a search with six unknown parameters: the gravitational wave amplitude 
ho, the polarization angle -0 (which depends on the position angle of the spin axis in 
the plane of the sky), the phase of the signal at a fiducial time 4>q, the inclination of 
the spin axis with respect to the line-of-sight i, the uncertainty in the absolute value 
of the signal frequency A/, and the frequency derivative A/. 

We use a Bayesian Markov chain Monte Carlo (MCMC) technique for this analysis 
as MCMC methods have been applied successfully to similar problems involving large 
numbers of parameters In a previous study [TI] . we used a Metropolis-Hastings 
(MH) algorithm ^|^] for a similar search, but with only five parameters (A/ being 
absent). When the frequency derivative A/ is included in the basic MH routine 
of ^31 the large correlation between A/ and A/ makes the parameter search difficult, 
and the basic MH algorithm becomes inefficient. In order to adequately sample 
the parameter space we implemented a combination of three different strategies for 
accelerating convergence of Markov chains: reparameterisation, the delayed rejection 
method of Tierney and Mira JS] (which is an adaptive version of the MH algorithm) , 
and simulated annealing 16J (which is a Monte Carlo approach to global optimization). 
The parameter Af is highly correlated with A/, and a strong correlation also exists 
between h$ and cos i. An initial transformation of these variables to near orthogonality 
yields a more tractable parameter space that is more effectively sampled. 

The heterodyne manipulation of the data used in this study is identical to that 
presented (by two of us) in an end-to-end robust Bayesian method of searching for 
periodic signals in gravitational wave interferometer data and is also described 
in 7 . A brief summary of this heterodyne technique is given in Sec. El Our delayed 
rejection method, as well as the reparameterisation strategy, is presented in Sec. |21 
In Sec. 0] we present the results of this study, using synthesized signals, for this 
six parameter problem. A brief discussion of the long term goals for this work are 
presented in Sec. [SJ 



2. The gravitational wave signal 

Gravitational waves from spinning neutron stars are expected to be weak at the Earth, 
therefore long integration periods are necessary to extract the signal. It is therefore 
important to take proper account of the antenna patterns of the detectors and the 
Doppler shift due to the motion of the Earth. 

As in previous studies (7J DZ| we consider the signal expected from a non- 
precessing triaxial neutron star. The gravitational wave signal from such an object is 
at twice its rotation frequency, / s = 2/ r , and we characterise the amplitudes of each 
polarization with overall strain factor, ho. The measured gravitational wave signal 
will also depend on the antenna patterns of the detector for the 'cross' and 'plus' 
polarisations, i*x,+, giving a signal 

s(t) = -F+(t;ijj)h a {l + cos 2 i) cos#(t) + Fx(*;V) ho cost sin \I>(i), (1) 
A simple slowdown model provides the phase evolution of the signal as 



where 



/ s (T-T ) + i/ s (T-T ) 



(2) 



T ' Tl 

T = t + St = t+ + AT. (3) 
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Here, T is the time of arrival of the signal at the solar system barycenter, (j> is the 
phase of the signal at a fiducial time T , f is the position of the detector with respect 
to the solar system barycenter, n is a unit vector in the direction of the neutron star, 
c is the speed of light, and AT contains the relativistic corrections to the arrival time 

m 

If f s and f s are known from (for example) radio observations, the signal can be 
heterodyned by multiplying the data by exp[— i\?(t)], low-pass filtered and resampled, 
so that the only time varying quantity remaining is the antenna pattern of the 
interferometer. We are left with a simple model with four unknown parameters h$, ip 1 
4>o and i. If there is an uncertainty in the frequency and frequency derivative then we 
have two additional parameters, the differences between the signal and heterodyne 
frequency and frequency derivatives, A/ and A/, giving a total of six unknown 
parameters. 

A detailed description of the heterodyning procedure is presented elsewhere 
I17j . Here we just provide a brief summary of this standard technique. The raw signal, 
s(t), is centered near twice the rotation frequency of the neutron star, but is Doppler 
modulated due to the motion of the Earth and the orbit of the neutron star if it is 
in a binary system. The modulation bandwidth is typically 10 4 times less than the 
detector bandwidth, so one can greatly reduce the effective data rate by extracting 
this band and shifting it to zero frequency. In its standard form the result is one 
binned data point, Bk, every minute, containing all the relevant information from 
the original time series but at only 2 x 1CP 6 the original data rate. If the phase 
evolution has been correctly accounted for at this heterodyning stage then the only 
time- varying component left in the signal will be the effect of the antenna pattern of 
the interferometer, as its geometry with respect to the neutron star varies with Earth 
rotation. Any small error, A/, in the heterodyne frequency will cause the signal to 
oscillate at A/ (plus the residual Dopper shift). We estimate the noise variance, u\, 
in the bin values, Bk, from the sample variance of the contributing data. It is assumed 
that the noise is stationary over the 60 s of data contributing to each bin. 



3. The adaptive Metropolis-Hastings algorithm 



After heterodyning, the signal on which we wish to carry out our MCMC analysis has 
the form [T7| 



V 



(t fc ;o) = ]F + (t k ;iP)h (l + cos 2 L)e iA *W - *-F x (t k ; ^)h cos^ A *« 



(4) 



where tk is the time of the fc^h bin Bk and a = (ho, ip, <fio, cos i, A/, A/) is a vector of 
our unknown parameters. A^(t) represents the residual phase evolution of the signal, 
equalling 4> + 2tt[A/(T - T ) + A/(T - T ) 2 /2]. The objective is to fit this model to 
the antenna output data 

Bk = y{tk\a) + e k , (5) 

where is assumed to be normally distributed noise with a mean of zero and known 
variance o\. Assuming exchangeability of the binned data points, Bk, the joint 
likelihood that these data d = {Bk} arise from a model with a certain parameter 
vector a is ^7] 



p(d\a) cx J^exp 

fc 



B k - V(t k ;a) 



o~k 



2 




[-X 2 (a)l 




= exp 


2 



(G) 
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where 



x 2 («) = E 



Bk - y(tk; a) 2 



(7) 



k 



In order to draw any inference on the unknown parameter vector a we need the 
(posterior) probability of a given d which can be obtained from the likelihood via an 
application of Bayes' theorem. The unnormalized posterior density 



is the product of the prior density of a, p(a), and the joint likelihood. Accordingly, 
appropriate priors have to be chosen for the particular parameters. In this study 
we use uniform priors with prior ranges [— 7r, 7r], [— 7r/4, 7r/4] and [—1,1] for the 
angle parameters <fio, if) and cost respectively. For ho we also specify a uniform 
prior with boundary [0,1000] in units of the rms noise |17|. For the frequency and 
spindown uncertainty we use suitable uniform priors with ranges of [— ^j] Hz and 
[-10~ 9 , 10~ 9 ] Hzs- 1 for A/ and A/, respectively. 

The normalized posterior density p(a\d) = p{a)p{d\a) / p(d) cannot be evaluated 
analytically, so we use Monte Carlo methods to explore p(a\d). If we can simulate 
from p(a\d), we can estimate all interesting quantities, including the posterior means 
of all parameters from the corresponding sample means, to any desired accuracy by 
increasing the sample size. 

However, drawing independent samples in a six-dimensional parameter space 
is not feasible. It has already been shown that MCMC methods can be used 
to parameterise gravitational wave signals of low signal-to-noise ratio |14j with 
four unknown parameters. These simulate a Markov chain, constructed so that 
its stationary distribution coincides with the posterior distribution and the sample 
path averages converge to the expectations. A minimal requirement for this is the 
irreducibility of the chain and hence the ability of the chain to reach all parts of 
the state space JT]. A specific MCMC technique is the MH algorithm H^] 
which does not require the normalization constant, only the unnormalized posterior 
density of Eq. (JSJ). We employed the MH algorithm for the four and five parameter 
pulsar detection problems ^3]. The efficiency of the MH algorithm depends heavily 
on the choice of the proposal density. Intuition suggests that the closer the proposal 
distribution is to the target, the faster convergence to stationarity is achieved. Default 
choices such as a Gaussian proposal or a random walk result in very slow mixing for this 
6-parameter problem. To increase the speed of convergence, we employed an adaptive 
technique, adaptive in the sense that it allows the choice of proposal distribution 
to depend upon information gained from the already sampled states as well as the 
proposed but rejected states. The idea behind the delayed rejection algorithm specified 
by 53 is that persistent rejection, perhaps in particular parts of the state space, 
may indicate that locally the proposal distribution is badly calibrated to the target. 
Therefore, the MH algorithm is modified so that on rejection, a second attempt to move 
is made with a proposal distribution that depends on the previously rejected state. 
This adaptive Monte-Carlo method JS] was generalized for the variable dimension 
case and renamed the 'delayed rejection method'. Since we have a fixed dimension 
problem here we implemented the original version P2] , and also the generalization [0?] 
that uses the reversible jump method. It turned out that the delayed rejection with 
the reversible jump method was not that beneficial for this particular problem and 
thus we will explain the original delayed rejection algorithm here. 



p{a\d) oc p(a)p(d\a) 



(8) 
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For the Metropolis-Hastings algorithm a new state in a Markov chain is chosen 
first by sampling a candidate a' from a certain proposal distribution g 1 (a'|a n ) usually 
depending on the current state a n and then accepting or rejecting it with a probability 
ai(a'\a n ) depending on the distribution of interest. This rejection is essential for the 
convergence of the chain to the intended target distribution. The choice of a good 
proposal distribution is important to avoid persistent rejections in order to achieve 
good convergence of a chain. However in different parts of the state space different 
proposals are required. When a proposed MH move is rejected, a second candidate a" 
can be sampled with a different proposal distribution q 2 (a"\a' , a n ) that can depend 
on the previously rejected proposal. Since a rejection suggests a bad fit of the first 
proposal, a different form of proposal can be advantageous in the second stage. To 
preserve reversibility of the Markov chain and thus to comply with the detailed balance 
condition, the acceptance probabilities for both the first and the second stage are given 
by EDI 



ai(a'\a n ) 
and 

a 2 (a"\a n ) 



min 1 



p{a')p(d\a')q l (a n \a') 
"' p(a n )p(d\a n )qi(a'\a n ) 



min 1, 



p{a")p{d\a") qi (a'\a") q2 {a n \a',a"){l- ai {a'\a")] 



(9) 



(10) 



p(a n )p(d\a n )q 1 (a'\a n )q 2 (a"\a n ,a')[l - ai(o'|o„) 

respectively. Fig. ^ illustrates the idea of delayed rejection. When the second stage 
proposal step is applied due to rejection of the first, the chain has, in order to 
preserve the reversibility, to imply a return path which comprises a fictive stationary 
Markov chain consisting of a fictive stage 1 proposal step from a" to a' which is 
rejected followed by an accepted fictive second stage move to a n ^H]- Although 



fictive 
Stage 2 
accepted 




Stage 2 
(timid) 



Figure 1. The delayed rejection method. In case of rejection of the first, bold 
step a second, more timid move is proposed. In order to maintain the reversibility 
of the Markov Chain the acceptance probability has to consider a fictive return 
path. 

the delayed rejection method provides better acceptance rates over the two stages, 
cross-correlations between the parameters still impede convergence of the Markov 
chain. Preliminary runs reveal that especially the parameters A/ and A/, and to a 
certain extent ho and cos i are highly correlated after the Markov chain has found a 
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potential mode. The consequence of which is poor mixing of the chain and therefore 
a reparameterisation is required. 

The coherence between A/ and A/ is obvious since the data is sampled from 
time tstart to iend; where the heterodyned signal traverses a frequency from /start = 
A/ + iA/-i star t to / ond = A/ + \ A/ ■ t cnd ; time t = is an epoch time when / = A/. 
Hence it is much more natural to work with these frequencies as new parameters and 
vary them with a certain correlation which influences A/ indirectly. The original 
parameters are then obtained by the simple linear transformation 

A/ = /start — ' ^ start 

and 

_ 2 /end ~ /start (12) 
^end ^start 

Since the Jacobian of this transformation is constant the prior distributions for the 
new parameters /start and / cn d are flat as well. 

Another cross-correlation can be observed between the parameters ho and cos l 
that arises from the fact that ho can be seen as a scaling factor and cos t as a nondinear 
weighting between the plus and cross polarisation part of the model. As seen in Eq. 
0J the plus part is multiplied by the factor a\ — i/io(l + cos 2 l) while the cross part 
encloses the term a 2 — 5/10 cos/,. The original parameters can be derived from 



h Q = 2 ^01 + - a\ \ , (13) 

and 

2«2 , v 

cost=— -. (14) 
ho 

As mentioned above, the prior distribution of the parameters ho and c = cos 1 are 
chosen uniform with joint probability density function 

f(ho,c) = \ *°^<fc«. -1<c<1, 

J v ' ; 10, otherwise, v ; 

where for this study lf l0 = 1000 in units of the rms noise. This implies a joint prior 
distribution for the parameters a\ and a 2 of the form 

m.n.a-,): j ■ if N<«i<^£^<%^ ]| detJ | (16) 




otherwise 



with Jacobian 



detJ=^^^=. (17) 

V a l ~ a 2 

Since the Jacobian is positive for the above restrictions we can write 

if \a 2 \ < ax < jr < 



9(ai,a 2 )={ hoV^' - ---- ii h0 -a' (18 ) 
0, otherwise. 

This joint prior density has the shape shown in Fig. [21 These reparameterisations result 
in a fast mixing Markov chain but still, the choice of a suitable proposal distribution 
is essential. Usually, a multivariate Normal distribution is utilized for the proposal 
distributions qi(a'\a n ) and q2(a"\a', a n ), with means equal to the current state and 
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a 2 



Figure 2. Joint prior density of a\ and for a given boundary Z/, for the 
parameter ho 

different variances depending on the stage. Larger variances are chosen for the 'bold' 
first stage steps, while smaller variances are more beneficial for the 'timid' second 
stage candidates. The covariance matrix has to comprise the correlation between the 
parameters /start and / on d since this correlation indirectly controls the parameter A/ 
as mentioned above. Hence choosing proposals for /start and / en d with a correlation 
of 1 would imply no change of A/ because both parameters are changed in the same 
way while a correlation of would have a great impact on A/ since /start and / cn d 
are changed completely uncorrclatcd. Thus the correlation between /start and / cn( j 
has to be treated randomly in order to control A/. Best results are obtained when 
a correlation of is chosen with probability 0.5 for the bold moves of A/ and a 
correlation of almost 1 otherwise for timid moves of A/. The proposals for the 
parameters a\ and a 2 are sampled independently since they represent scaling factors 
for the plus and cross polarisation part, respectively. Finally we have to consider the 
correlation between the original parameters ip and <f>o which are not reparameterised. 
Pilot runs show that they are highly correlated. Hence the proposal distribution is 
adapted accordingly. 

Unfortunately, the posterior distribution features very narrow modes in a large 
parameter space that has to be scanned. Thus a simple Normal distribution is not 
suitable for a proposal distribution as pilot runs have revealed. Instead, a proposal 
distribution with long tails and strong narrow mode is required. This can easily be 
achieved by generating a random sample between two boundaries bi and bh for the 
standard deviation of the proposal by generating a random weight for the weighted 
geometric mean of these two boundaries. Hence we sample standard deviations 
according to a = b^b\~ w where w ~ (3(a, b) is Beta-distributed with parameters a 
and b. The resulting proposal distribution is symmetric with very long tails and 
a strong narrow mode. In order to obtain higher standard deviations for the first 
stage the choice of w ~ (3(2, 1) (with mean |) is adequate while for the second stage 
w <~ /3(1,2) (with mean 1/3) samples smaller standard deviations. 

The implementation of the ideas outlined above leads to reasonable acceptance 
rates and hence to a much better convergence of the Markov chain. While during the 
burn-in period it is mainly the stage 1 candidates that are accepted, the Markov chain 
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is driven mainly by stage 2 candidates after the burn-in. But still, the stationary 
distribution features many distinct modes that carry the risk of trapping the Markov 
chain. Therefore, we regard the posterior as a canonical distribution 



p{a\d) oc p(a)p(d\a) ocp(a)exp 

X 2 (a)-21ogb(a)] 



oc exp 



oc exp [-P (x 2 (a) - 2 log [p(a) 



(19) 



with inverse temperature (3. During the burn-in period this inverse temperature can 
pass through values starting at a low value (thus high temperature) and ending up 
at (3 = | which coincides with the posterior distribution. This simulated annealing 
technique was introduced by Metropolis et al. ^5] an d allows scanning of the whole 
parameter space by permitting larger steps. For the annealing schedule an exponential 
temperature curve is applied. For a certain number of iterations t s , it starts with an 
inverse temperature (3o until it reaches (3 = h. The inverse temperature follows the 
function 



I3(t) 



/3 exp flog 



if < t < t s 
if t > 



(20) 



depending on the current iteration t. Since the starting temperature is dependent on 
the data set which is influenced by the amplitude kg of the signal it has to be adapted 
accordingly. 



4. Results with simulated signals 

We have synthesized fictitious data, and passed it through our six parameter MCMC 
routine. The presentation of results here is similar to that of the four and five 
parameter study of |14|. The artificial signals were embedded within white and 
normally distributed noise. The ability of the MCMC algorithm to successfully find 
the signal and estimate the six parameters was demonstrated, and is presented below. 
The artificial signals s(t) were synthesized assuming a source at RA = 4 h 41 m 54 s and 
dec = 18° 23' 32", as would be seen by the LIGO-Hanford interferometer. The signals 
were then added to noise; we assumed a signal at 300 Hz and a corresponding noise 
spectral density of at that frequency of h(f ) = 8 x 10 _23 Hz -1 / 2 . The amplitude of 
the signal used in our test runs was varied in the range ho = 4.0 x 10~ 24 to 1.5 x 10 -22 . 
The length of the data set corresponded to 14400 samples or 10 days of data at a rate 
of one sample per minute (which was the rate used for the LIGO/GEO SI analysis 
described in [7]). 

In Fig. |3 we display the MCMC generated posterior probability distribution 
functions (pdfs) for an example signal. The real parameters for this signal were: 
h = 1.5xl0~ 22 , V = 0.4, 0o = 1.0 (both in radians), cost = 0.878, A/ = 7.0xlO~ 3 Hz 
and A/ = —2.5 x 10~ 10 Hzs -1 . For this example the program ran for 10 6 iterations. 
For a signal this large only about 2.5 x 10 4 iterations were needed for the burn-in, 
and this data is discarded from the analysis. Short-term correlations in the chain were 
eliminated by 'thinning' the remaining terms; we kept every 250 th item in the chain. In 
this example the MCMC yielded median values and 95% posterior probability intervals 
of h = 14.91 x 10~ 23 (13.41 x 10~ 23 to 15.84 x 10~ 23 ), V = 0.439 (-0.552 to 0.707), 
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cf> = 0.964 (0.696 to 1.958), cosi = 0.884 (0.828 to 0.988), Af = 6.99999772 x 10" 3 Hz 
(6.99999217 x 10~ 3 Hz to 7.00000314 x 10~ 3 Hz) and Af = -2.4999541 x 10~ 10 Hz s" 1 
(-2.5000767 x lO^Hzs- 1 to -2.4998272 x lO-^Hzs" 1 ). The 95% posterior 
probability interval is specified by the 2.5% and 97.5% quartiles of £>(a,|d). 





Figure 3. MCMC estimates of the posterior pdf (kernel density) for the 
six parameters ho, ip, <j>0, cosl, Af and Af. This synthesized signal had 
real parameters of: ho = 15 X 10 -23 , tp = 0.4, (f>o = 1.0, cosl = 0.878, 
Af = 7.0 X 10~ [i Hz and Af = -2.5 X 10 -10 . The mean of the ho distribution 
here is 14.84 X 10 -23 . 

With the noise level used, h(f) — 8 x 10~ 23 , we were able to successfully detect 
signals with amplitudes of ho > 4.0 x 10~ 24 with 10 days of data. This should be 
compared with the results presented in |14| where with just four parameters (ho, 
if), 4>o and cosl), we were able to confidently detect signals with an amplitude four 
times smaller. The addition of the new frequency parameters has the disadvantage of 
complicating the search due to the corresponding increase in the size of the parameter 
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space. For our study, we let the initial burn-in of the Markov chain last for as long 
as 3.5 x 10 5 iterations, and if the signal was not found by this time the search was 
terminated. It may be possible to find smaller signals with a longer burn-in. 

Fig. 0| shows the MCMC estimated posterior for the smallest value of the 
parameter ho that we were able to identify with the MCMC code. The true parameter 
values for this run were ho = 4.0 x 10~ 24 , tp = 0.4, </> = 1.0, i = 0.5 (cost = 0.878), 
A/ = 7.0 x 10~ 3 Hz and A/ = -2.5 x 10 -10 Hzs -1 . In this run the MCMC 
yielded a mean value and 95% posterior probability interval of ho = 4.8 x 10 -24 
(0.34 x 10~ 24 to 0.74 x 10~ 24 ). Fig. displays the MCMC estimated posterior 
for the parameters A/ and A/, which provides mean values and 95% posterior 
probability intervals of A/ = 7.0 x 10~ 3 Hz (6.9998 x 10~ 3 Hz to 7.0002 x 10~ 3 Hz), 
and Af = -2.500 x 10- 10 Hz,s _1 (-2.505 x 10" 10 Hz.s" 1 to -2.496 x 10- 10 Hz,s _1 ). 
As can be seen from Figs. 0] and 03 even with small signal level it is still possible to 
extract the most astrophysically important parameters. For this MCMC run there 
were a total of 10 6 iterations, with the first 3.5 x 10 5 as the burn-in. 



Figure 4. MCMC estimate of the posterior pdf (kernel density) for the parameter 
ho from a six parameter search using synthesized data. The real parameter value 
for this signal was ho = 4.0 x 1CT 24 . This was the smallest signal detectable by 
the MCMC method for the noise level used. 



5. Discussion and conclusions 

In the simplest application, the method demonstrated here could complement searches 
for signals from known pulsars I17| : our method could be used to verify the 
frequency and frequency derivative values. The real advantage of the technique would 
come about in a search for a signal at a known location, but where the frequency 
information pertaining to the neutron star is not well known; a search for a signal 
from SN1987A ^Uj would be a possible application. In the demonstration here the 
heterodyning process provides a band of 1/60 Hz. It would be straightforward to 
expand this search to a bandwidth of 5 Hz by running the code on 300 processors, a 
task easily accomplished on a cluster of computers. For 10 days of data it takes a 
single 2.8 GHz personal computer approximately an hour to conduct about 3.3 x 10 4 
iterations of our MCMC code. There are more iterations done per time interval at the 
beginning of a run because at that time more stage-1 steps are accepted. We believe 
that these MCMC methods offer great potential benefits for gravitational radiation 
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0.0069996 0.0069998 0.0070000 0.0070002 0.0070004 -2.510e-10 -2.505e-10 -2.500e-10 -2.495e-10 



N = 2601 Bandwidth = 1 .933e-08 N = 2601 Bandwidth = 4.258-1 4 

Figure 5. MCMC estimate of the posterior pdfs (kernel densities) for the 
parameters Af and Af from a six parameter search using synthesized data with 
the smallest detectable signal ho = 4.0 X 10 — 24 . The real parameters for this 
signal were: Af = 7.0 X 10~ 3 Hz and Af = 2.5 X 10" 10 Hz.s" 1 . 

searches where the signals depend on a large number of parameters. 
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